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THE  ITERATIVE  DETERMINATION  OF  MODEL 
PARAMETERS  BY  NEWTON'S  METHOD 


The  problem  that  we  would  like  to  solve  is  to  determine  numbers 
aa,  aa,  ....  a„  and  An  A2,  ...»  An  so  that,  with  regard  to  an  ob- 
servation e{t)  and  for  unknown  N , the  equation 

N 

e(t)  a*  s(t-^K  ) (1) 

K = 1 


is  satisfied. 

Let  us  define 

N 

y(t;x)  x*  s(t-xN+K  ) - e(t).  (2) 

K =1 

In  keeping  with  the  spirit  of  Newton's  method  for  finding  the  root  of 
a function  of  one  variable,  we  are  looking  for  a vectorial  increment 
fix  such  that 


S N 

y(t;x+fix)  = y(t;x)  +V  ^ fix*  + terms  of  higher  order  =0. 

i-J  0Xk 
K =1 

Thus  we  have  an  equation 
a n 

S 6Xk  = (3) 

K = X 

from  which  to  determine  the  incremental  vector  components  fix*  . 

At  thi3  point  it  is  important  to  recall  that  both  members  of  Equation  3 
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are  functions  of  t.  In  order  to  eliminate  t and  obtain  2N  equation 
in  2N  unknowns  we  could,  for  example,  choose  2N  values  of  t, 
t,,  t2,  ....  t?  N and  evaluate  both  members  of  Equation  3 at  these 
t values.  In  a real  situation  this  procedure  has  the  drawback  that 
the  demands  for  accuracy  placed  on  the  values  of  y and  its  partial 
derivatives  at  the  special  t points  would  be  too  great.  In  the  presence 
of  noise  (due  to  inaccuracies  in  measurement,  in  the  model,  round-off 
errors,  etc.  ) the  above  method  of  creating  2N  equations  in  the  2N 
unknowns  should  be  rejected  as  being  unstable. 

Actually,  however,  this  method  is  just  one  realization  of  a general 
method  of  deriving  2N  equations  in  2N  unknowns  from  Equation  3. 

The  general  method  postulates  that  2N  linear  functionals  be  chosen 
y i,  y 2,  ...,  y2N  and  applied  to  both  members  of  Equation  3.  The 
result  is  then 

S N 

S Vl(I?;)6x'  =-Vi-(y)-  (L  =1'2 2N)*  (4) 

K = 1 

We  now  have  a non-homogeneous  linear  system  of  equations  in  the 
unknown  6xx  (k  = 1,  2,  ...,  2N)  and  this  system  has  a unique  solution 
only  if  the  matrix  M with  elements 


is  non- singular  and  the  vector  g with  elements 

-yL  (y)  (6) 

does  not  vanish. 

It  is  necessary  to  say  a word  about  the  notion  of  a linear  functional 
y . y associates  with  each  member  f of  a class  of  functions  a 
number  y(f)  in  such  a way  that  if  f and  g are  any  two  functions  in 
the  class,  a and  0 are  any  two  numbers,  and  that  if  the  function 
0(1  + fig  also  belongs  to  the  class,  then 


y (<*£  + 0g)  = a y (*)  + 0 y (g). 
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If  the  class  of  functions  consists  of  all  those  which  are  square  in- 
tegrable  over  the  real  axis,  then  by  the  celebrated  Riesz-Fischer 
theorem*  every  bounded  functional  y has  the  representation 

00 

y (0  = j i(x)  h (x)  dx  , 

-CO  ^ 

where  h is  a square  integrable  function  corresponding  to  the  linear 
functional  y and  the  bar  denotes  complex  conjugate. 

After  this  brief  digression  let  us  return  to  Equation  4 to  see  how 
it  can  be  used  to  approach  a solution  by  iteration.  Suppose  at  the  *TH 
step  an  approximate  solution  ZR  has  been  calculated.  To  determine 
the  next,  hopefully  improved,  approximant  ZR,  we  set 


6x*  = Z*+1 

-ZR 

and  let 

0 

Ml/  s-yL" 

( *Z\  and 
V&x,  ) 

(7) 

8lr  s-y/ 

(y*)  . 

(8) 

Then 

M"  ZR+1  = 

MrZ"  -gR 

(9) 

and  if  M*  has  an  inverse,  then 

Z*+l  = Z*  - (M")  "l  gR.  (10) 

It  can  be  shown-*-  that  if  the  initial  guess  Z°  is  sufficiently  close 
to  a solution,  then  the  sequence  of  approximants  Z°,Z1,Za,  ...  con- 
verges to  the  solution. 


* Halmos,  Paul  P.,  Introduction  to  Hilbert  Space,  Chelsea 
Publishing  Company,  New  York,  1951,  p 31. 

t Stern,  M.  L. , "Sufficient  conditions  for  the  convergence  of  Newton's 
Method  in  complex  Banach  spaces",  Proc.  Amer,  Math.  Soc.,  vol  3, 
(1952)  pp  858-863.  “ 


67TMP-64 


5 


& 


It  is  apparent  that  one  of  the  major  problems  with  the  method  is 
the  selection  of  a good  initial  approximant  Z°.  Two  methods  have  been 
been  devised  to  obtain  approximate  solutions  for  Equation  1.  The 
first  will  be  called  the  Fourier  method.  Taking  the  Fourier  transform 
of  both  members  of  Equation  1,  the  result  is 


E(f) 


2 


aK 


- 2 Trif  An 

0 


S(f). 


(11) 


This  expression  i6,  however,  valid  only  for  those  frequencies  for 
which  S(f)  and  E(f)  do  net  vanish,  that  is,  only  for  the  frequencies 
that  lie  in  the  common  band  in  case  s and  consequently  e are 
bandlimited.  If  in  Equation  11  we  divide  both  members  by  S(f)  and 
then  take  the  inverse  Fourier  transform  of  the  result,  we  obtain 


d(t) 


N 


n 


(12) 


where  ft  represents  the  pertinent  band  of  frequencies. 


Ifft=|f-^-sf£~-  , then 


d(t) 


sin  ffw(t-AK  ) 

rr(t-AK ) 


(13) 


and  the  values  of  Ax,  As.  ....  An  can  be  estimated  as  the  position 
of  the  maxima  of  d(t).  The  values  of  ax,  a2,  . . .,  aN  are  approx- 
imated by  the  magnitudes  of  d(t)  at  the  observed  maxima.  The  limit 
of  resolution  of  this  method  is  about  1/w,  i.  e.,  if  A.  - A-  < 1/w, 

1 lTl 

then  the  corresponding  two  maxima  will  have  moved  together  so  as  to 
yield  only  one  maximum. 


< f < - wc 


< f < wc 


I 


\ 
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} 


i 


ft 


i. 


then 


d(t)  -^T  aK 


sin  it  (2we  +w/2)(t-A,c  )-sin  v (2we"w/2)(t-A|(  ) 


(14) 


K =1 


or 


d(t)  = aK  cos  2irwe(t-AK)  Sl;;r^2)(t-^  * 


(15) 


K =1 


w w 

— r £ f £ W.  + — 

2 c 2 


Again  the  same  method  for  estimating  the  A's  and  a's  can  be  used. 

If  n = | f | wc  - 

then 

affiwc(t-AK)  sinffw(t-An ) 


d(t) 


a,  e 


ff(t-AK ) 


(16) 


* =i 


i 

i 


and  the  real  part  will  again  have  maxima  near  the  points  t-A  and  the 
magnitudes  of  these  maxima  will  be  close  to  a* . 

The  second  method,  based  on  correlation,  can  also  be  used  for  obtain- 
ing first  estimates.  Let  us  multiply  both  members  of  Equation  1 by 
s(t-z)  and  integrate  over  all  values  of  the  variable  t . We  can  then  write 

N 

Rs.e  (z)  = £ an  Rs,S  <*-A,c),  (17) 

* =1 

where 

CD 

Rf  ,a  <*)  = f Kt-z)  g(t)dt  . (18) 

J 

• 00 

If  the  Ak  are  sufficiently  separated,  i.  e.,  by  more  than  i/w,  then 
the  position  of  the  maxima  of  R;  ,E  (z)  provides  first  estimates  for 
Ax , As  > . . . An  and  the  magnitudes  of  the  maxima  divided  by  Rs  t , (o) 
provide  estimates  for  at,  ag,  a„. 
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Having  found  first  estimates  of  , Aa»  ...»  AN  and  aj.  aa»  . .., 
aN,  we  must  now  seek  to  improve  them  by  iteration.  The  requirements 
for  a good  iterative  method  are  : 

1.  Improvement  of  resolution.  If  one  maximum  actually 
corresponds  to  two  separate  components,  the  method  should 
be  able  to  resolve  the  two  components. 

2.  If  some  secondary  maxima  have  been  confounded  with 
primary  maxima  the  iterative  method  should  be  capable  of 
eliminating  them. 

3.  Finally,  an  improvement  in  accuracy  in  the  evaluation 
of  Ai»  ...,  An  should  result  from  the  application  of  the 
method. 

Once  the  initial  estimates  have  been  made  we  turn  to  Newton's 
method  for  their  improvement.  To  this  end  we  must  choose  functionals 

y i • * • • » ys n • 


Let 


yi(f)  = l f(t)  C(£-Xi+N) 


dt 


and 


U)  = j ^ l-^a,(t-X.+N)}dt  , 


where  s'  denotes  the  derivative  of  s.  The  matrix  M",  defined  by 
Equation  7,  takes  the  form 


00 

= J s(t-xJ+N)  s(t-xL  +l)  dt 

••CD 


M 


« 

l+N* 


8<t-x5+N>  ^ "XJ  s'(t-xL+*))  dt 


-yv 
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{.-xj  s'(t-x"+n)l  s(t-x^+N)dt 


- / i-*;  •■(*-*:«>  > s'(t'*:+J  1 d* 

~co 

for  L * U,  ...»  N and  K = 1,  2,  ...»  N- 

I,  • us,  „ notioi.  of  convolution  or  correlation  we  can  write  the 
above  n.  ..  x\  * .emente  more  conveniently  as 


m;„  = 

Rs , S (XK  +N 

Xt  +N  ^ 

4-n,  k 

= -X.:  Rs<5. 

<X*+N 

- X?+N) 

^ U , X +N 

= -XI  Rs',s 

(A  +N 

- xj  fN) 

MJ+n<k+n  = x2xi  R»*,s'  (xk+h-x"+n) 

for  L = 1.2,  ...»  N and  K = 1,  2 N* 

But 


dR,,«(T)  _ 

b T 


f s'(t-r)  S(t)  dt  = -R s' , s It) 


^R5,»  tr) 

»T* 


= | s'(t-r)  » '(t)  dt 


Rj'  , s' 
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Therefore, 


M[iK  = R..,  (T) 


T = (XK  +N  "Xl  + N ) 


ml*+n,k  = x* 


mr  x - x* 

1V1L  *K+N  aK  ^ 


T = (*I  +N  - x*  +N  ) 


T = -(xJ+K  -xf+M) 


T = (XK+N  ”XL+n) 


(19) 


for  L = 1,2,  ...,  N K = 1,2,  ....  N K = 1, 2,  . . . , N 

The  vector  g*  whose  components  are  given  by  Equation  8 takes 
the  form 


00 

gj  = J e(t)  s(t-x*+N)  dt  = Rs,e  (x£+n) 


00 

g?+N  = J e(t)  { -x^  s'(t-x*+H)  } dt  = xL  -^.E 
-00 


(20) 


r=x!!+ 


Equations  19  and  20  show  that  the  only  data  that  we  require  for  the 
application  of  Newton's  method  are  the  correlation  functions  Rs,s 
and  Rs>e.  If  these  are  furnished  by  observation,  numerical 
differentiation  will  yield  the  remaining  matrix  and  vector  elements. 

In  the  course  of  iteration  it  is  necessary  to  update  the  matrix 
M and  the  vector  g at  each  step.  This,  however,  is  very  easily 
done  since  it  involves  no  more  than  a table  look  up  of  the  values  of 
functions  and  their  derivatives  at  new  values  of  their  arguments. 


To  test  the  ideas  presented  above  by  means  of  an  example,  the 
choice  was  made  of 


S(t) 


8 in  (fft) 
(irt)  * 
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This  choice  is  motivated  not  only  by  the  desire  to  generate  a test, 
but  also  by  the  notion  that  upon  transformation  by  the  Fourier  method 
Equation  1 takes  the  form  of  Equation  13  which  is  the  same  as  that 
of  Equation  1 when  the  above  choice  is  made.  Therefore,  with  this 
choice  we  can  proceed  to  solve  Equation  13  directly  by  Newton's 
method. 

To  determine  the  matrix  elements  consider 


R5.s<t)  = J 


ain  n (t-r)  Sin  n t 
it  (t-r)  fft 


00 

J es*i(T 

-CO 

X(f) 

2 df  sin  nr 


ITT 


where  X(f)  = 1 for  - 1/2  £ f s 1/2  and  X(f)  = 0 for  |f|>  1/2. 


dRs . s xcosx  - smx 

— - — < — - 

3t  xz  x = it  r 


-2xa  cos  x + x(2-xa)  sin  x 

A program  (listed  in  the  Appendix)  was  written  to  carry  the  cal- 
culations through  numerically.  Tests  were  performed  on  a weighted 
sum  of  delayed  and  truncated  replicas  of  six  irt/irt.  Then  tests  showed 
that  convergence,  when  it  occurred,  was  rapid.  We  had  numerical  dif- 
ficv’lties  in  high  dimensions  because  of  the  many  matrix  inversions 
invclved.  The  program  was,  however,  al  e to  resolve  two  pulses  that 
were  separated  by  0.  2,  that  is,  by  0.  2 of  the  normal  resolution  limit. 

In  order  to  facilitate  decisions  as  to  the  best  value  of  N,  the  mean 
square  difference  between  the  function  e and  its  approximant  was 
calculated  and  that  solution  was  adopted  as  the  final  one  which  had  a 
number  of  components  yielding  the  smallest  mean  square  difference. 

A possible  improvement  might  be  had  by  evaluating  the  gradient 
of  y,  i.e.,  { By/bx^  },  at  some  point  other  than  the  last  approximant 
to  the  zero.  Work  to  realize  this  idea  is  going  on  at  present  and  will 
be  reported  subsequently. 


aaRs,»  _ [ 

dr3  x4  l 
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CONCLUSION 

A method  for  solving  the  functional  equation 
a 

e(t)  = ^ ax  s{t-AK  ) 

*=i 

is  presented.  The  method  adapts  the  idea  of  Newton's  method  to  the 
case  at  hand  which  differs  from  the  classical  situation  in  that  the 
functions  whoso  zeroes  are  sought  are  themselves  elements  of  a 
function  space.  To  illustrate  the  method  a program  was  written 
that  demonstrated  a marked  improvement  in  resolution  for  pulses 
of  the  form  sin(trt)  / (fft).  One  of  the  drawbacks  of  the  methoc  is  that 
initial  guesses  have  to  be  close  to  the  actual  solution.  Another 
difficulty  arises  from  the  fact  that  since  N is  unknown,  solution 
must  be  attempted  with  various  values  of  N. 


-"*** ...w-wrtiN, 
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APPENDIX 


0«18  SB  MON  09/0S/67 


1 FILE  K-R 
10  DIM  E<256) 

20  DIM  M(20»  20)» 9(20* 20 >* F (20* I )»Q(20» 1)*A(10)*D(10) 

25  READ  TO* AO* BO* B 

30  FOR  1*1  TO  TOSREAD  FILE  1*E(I>\NSXT  I 
40  Pi=3* 1415927 
50  FOR  S*1  TO  100\ INPUT  T1.L2 
70  T=24T1 

110  FOR  I->1  TO  T1MNPUT  A(i  >*  DU  >\NEXT  I 
120  FOR  LI  M TO  L2 
130  MAT  M«Z£R<T*T> 

140  MAT  V*Z£R(T*T> 

150  MAT  Q=ZER(T* 1 > 

200  MAT  F=Z£R(T* 1 ) 

210  FOR  U*1  TO  TtVFOR  1=1  TO  TO 
220  A1=P1*<A0*1-B0“D( J) ) 

230  IF  ABS(A1 ><«. OOOOOt  THEN  270 
240  F<J* 1 >=FC J» 1 ) + E( 1 )*SIN(A1 )/ <A1*B) 

250  F<T1+J*  1 >=F(T1«-J*  1 >-E(  I >*P1*<A1*C0S(A1  )-SlN(Al  ) )/  (A1  «2» 

2f0  00  TO  280 

270  F(J»  1 >«F( J*  1 >♦£<  I )/B 

280  NEXT  1 

285  F(T1*J* I )=A(J)*F(T1*J» 1 ) 

290  NEXT  J 
295  MAT  F--i«0>*F 
300  FOR  1=1  TO  T1 
310  FOR  J=1  TO  T1 
320  K«P1*<D(J>-D(I>) 

325  IF  ABS<KX=.0000l  THEN  400 
330  M<l»J)=SlN<K)/(k*B) 

340  H<T1M»T1**J1=A< 1 )*A< J)*<P1 *?>*B*(2=K*C0S(K)M  CK»2>-8>*SIN(K>)/<K»3) 
350  M( l*Tl+J)=Pi *A(U)* (K*COS<K)“SIN<K>  >/ (K»2  > 

380  GO  TO  490 
400  M<l*J)«i/& 

410  M<T1* l»Tl*J)«A( I)*A(J>*B*(Plt8)/3 
420  M<1.T1*U>=0 
430  M<THI*J>=0 
490  NEXT  J\NEXT  1 

495  FOR  1=1  TO  Tl'.FOR  J=1  TO  Tl\M<I*Tl* J>=M(J*T»*I> 

496  NEXT  J\NEXT  I 

497  FOR  1=1  TO  T\FCR  J=l  TO  TMF  ABS(M(  1* J> )>=  10  EMI  THEN  500 

498  M(I*J>=OVNEXT  JSNEXT  I 
500  MAT  W=INV<M) 

510  MAT  Q=U*F 
520  FOR  0=1  TO  T1 
530  A<U)»Q(U*1> 

540  D(J)=D(J>+Q(T1+J* I > 

SSO  NEXT  J 
S60  NEXT  LI 
600  FOR  1=1  TO  Tt 
610  PRINT  UAdHDCt) 

620  NEXT  1 
630  NEXT  S 
700  DATA  100* <1*5.1 
9999  END 
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